library(flowCore)
library(flowStats)
library(flowViz)
library(ggplot2)
library(lessR)
##install.packages("lessR")
library(gridExtra)
##library(multicore)
source("/home/haseong/alpha/dev/FACS/facs_analysis_library.R")

if(.Platform$OS.type=="windows"){
  root<-"D:\\mydocs\\2013\\dev\\ProteinSW\\"
  data.dir<-"D:\\mydocs\\2013\\dev\\ProteinSW\\"
}else{
  root<-"/home/haseong/alpha/dev/FACS/IspS/"
  data.dir<-"/home/haseong/alpha/dev/FACS/DmpR/20131005_dmpR_E01_400_700/"  
  data.dir2<-"/home/haseong/alpha/dev/FACS/DmpR/20131008_dmpR_E01_400_700/"  
}
setwd(root)



x<-list()
cnt<-1

##=========
data.dir<-paste(root, "110826", sep="")
##file.names<-dir(data.dir, pattern="[0-9]$")
file.names<-dir(data.dir, pattern="GESS") ## need to check the dataset name starts with "GESS"
## read info
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
x<-list()
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-tmp.info
##=========
data.dir<-paste(root, "110828", sep="")
file.names<-dir(data.dir, pattern="GESS")
## read info
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
##data.info<-tmp.info
data.info<-rbind(data.info, tmp.info)
##=========
data.dir<-paste(root, "110830", sep="")
file.names<-dir(data.dir, pattern="GESS")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
##data.info<-tmp.info
data.info<-rbind(data.info, tmp.info)
##=========
data.dir<-paste(root, "110920", sep="")
file.names<-dir(data.dir, pattern="IspGESS")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
##data.info<-tmp.info
data.info<-rbind(data.info, tmp.info)
##=========
data.dir<-paste(root, "110927", sep="")
file.names<-dir(data.dir, pattern="IspGESS")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
##data.info<-tmp.info
data.info<-rbind(data.info, tmp.info)

##=========
data.dir<-paste(root, "110929", sep="")
file.names<-dir(data.dir, pattern="IspGESS")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
##data.info<-tmp.info
data.info<-rbind(data.info, tmp.info)

##========= various nitrotoluene derivatives
data.dir<-paste(root, "111107", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info2.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

##========= Sensitivity comparison of Toluene, Benzene, Isoprene 
data.dir<-paste(root, "111123", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

##========= Sensitivity comparison of Isoprene 
data.dir<-paste(root, "111123", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

## ============== 120228 various cell strain
data.dir<-paste(root, "120228 various cell strain", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

## ============== 121204 IspGESS promoter1
data.dir<-paste(root, "121203 IspGESS promoter1", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

## ============== 121204 IspGESS promoter2
data.dir<-paste(root, "121204 IspGESS promoter2", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

## ============== 130122 B.subtilis candidates
data.dir<-paste(root, "130122 B.subtilis candidates", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)

## ============== 20130307 Bindu sensitivity
data.dir<-paste(root, "20130307", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info[,-24])
data.info[which(data.info=="")]<-NA


## --- plot

tick_locations = c(0, 187.9363, 256, 434.9363, 512, 768, 1024)  
tick_labels = round(10^(tick_locations/256))


dh5a.idx<-which(data.info[,4]!="BL21(DE3)")
bl21.idx<-which(data.info[,4]=="BL21(DE3)")

tmp<-facs.summary.plot(x[[5]], ncells=10000, view.plot=T)
savePlot("DH5a_example.jpg", type="jpeg")
tmp<-facs.summary.plot(x[[146]], ncells=10000, view.plot=T)
savePlot("bl21_example.jpg", type="jpeg")


layout(matrix(c(1:4), 2, 2, byrow = FALSE))


## fsc dh5a
fsc.arr<-c()
cnt<-1
for(i in dh5a.idx){
  tmp<-facs.summary.plot(x[[i]], ncells=10000, view.plot=F)
  fsc.arr<-c(fsc.arr, tmp$sel.xx[,"FSC-H"])
  cat("---", cnt, "/", length(dh5a.idx), "---\n")  
  cnt<-cnt+1
}

##plot(density(10^(ssc.arr/256), width=20))
plot(fscd<-density(fsc.arr, width=10), xlab="FSC", xaxt = "n", xlim=c(0,1024), bty="n")
axis(1, at = tick_locations, labels = tick_labels)
fsc.max<-fscd$x[which.max(fscd$y)]
abline(v=fsc.max, lty=2)
text(fsc.max, max(fscd$y)*0.95, round(10^(fsc.max/256),2), cex=2, pos=4)
fsc.bottom<-10^(quantile(fsc.arr, 0.05)/256)
fsc.up<-10^(quantile(fsc.arr, 0.95)/256)


## ssc dh5a
ssc.arr<-c()
cnt<-1
for(i in dh5a.idx){
  tmp<-facs.summary.plot(x[[i]], ncells=10000, view.plot=F)
  ssc.arr<-c(ssc.arr, tmp$sel.xx[,"SSC-H"])
  cat("---", cnt, "/", length(dh5a.idx), "---\n")  
  cnt<-cnt+1
}

##plot(density(10^(ssc.arr/256), width=20))
plot(sscd<-density(ssc.arr, width=10), xlab="SSC", xaxt = "n", xlim=c(0,1024), bty="n")
axis(1, at = tick_locations, labels = tick_labels)
ssc.max<-sscd$x[which.max(sscd$y)]
abline(v=ssc.max, lty=2)
text(ssc.max, max(sscd$y)*0.95, round(10^(ssc.max/256),2), cex=2, pos=4)
ssc.bottom<-10^(quantile(ssc.arr, 0.05)/256)
ssc.up<-10^(quantile(ssc.arr, 0.95)/256)



##================================================


## fsc dh5a
fsc.arr<-c()
cnt<-1
for(i in bl21.idx){
  tmp<-facs.summary.plot(x[[i]], ncells=10000, view.plot=F)
  fsc.arr<-c(fsc.arr, tmp$sel.xx[,"FSC-H"])
  cat("---", cnt, "/", length(bl21.idx), "---\n")  
  cnt<-cnt+1
}

##plot(density(10^(ssc.arr/256), width=20))
plot(fscd<-density(fsc.arr, width=20), xlab="FSC", xaxt = "n", xlim=c(0,1024), bty="n")
axis(1, at = tick_locations, labels = tick_labels)
fsc.max<-fscd$x[which.max(fscd$y)]
abline(v=fsc.max, lty=2)
text(fsc.max, max(fscd$y)*0.95, round(10^(fsc.max/256),2), cex=2, pos=4)
fsc.bottom<-10^(quantile(fsc.arr, 0.05)/256)
fsc.up<-10^(quantile(fsc.arr, 0.95)/256)

## ssc bl21
ssc.arr<-c()
cnt<-1
for(i in bl21.idx){
  tmp<-facs.summary.plot(x[[i]], ncells=10000, view.plot=F)
  ssc.arr<-c(ssc.arr, tmp$sel.xx[,"SSC-H"])
  cat("---", cnt, "/", length(bl21.idx), "---\n")  
  cnt<-cnt+1
}

##plot(density(10^(ssc.arr/256), width=20))
plot(sscd<-density(ssc.arr, width=20), xlab="SSC", xaxt = "n", xlim=c(0,1024), bty="n")
axis(1, at = tick_locations, labels = tick_labels)
ssc.max<-sscd$x[which.max(sscd$y)]
abline(v=ssc.max, lty=2)
text(ssc.max, max(sscd$y)*0.95, round(10^(ssc.max/256),2), cex=2, pos=4)
ssc.bottom<-10^(quantile(ssc.arr, 0.05)/256)
ssc.up<-10^(quantile(ssc.arr, 0.95)/256)



savePlot("ssc_fsc.jpg", type="jpeg")


























myfplot<-function(wf, data.info, ctl.idx, ncells=10000, xlimi=c(0,800)){
  n.exp<-length(wf)
  par(ps=22, mfrow=c(n.exp,2), mai=c(0.4, 0.4, 0.3, 0.2), oma=c(0,0,0,0))
  ## control 
  xx<-exprs(wf[[ctl.idx]])
  del.idx<-which(xx[,"SSC-H"]<=50)
  xx<-xx[-del.idx,]
  ## select 10000 cells from
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  sel.idx<-o[1:ncells]
  sel.xx<-xx[sel.idx,]
  cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), i, "\n");flush.console()
  plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  ##plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", main=data.info[ctl.idx,23],xlim=c(0,500))
  plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", xlim=xlimi, bty="n", main="")
  text(0, mean(density(sel.xx[,"FL1-H"])$y), paste(data.info[1, "effector"], data.info[1, "effector.level"], sep=" "), pos=4)
  for(i in c(1:n.exp)[-ctl.idx]){
    xx<-exprs(wf[[i]])
    del.idx<-which(xx[,"SSC-H"]<=50)
    xx<-xx[-del.idx,]
    ## select 10000 cells from
    fsc.density<-density(xx[,"FSC-H"])
    fsc.max<-fsc.density$x[which.max(fsc.density$y)]
    ssc.density<-density(xx[,"SSC-H"])
    ssc.max<-ssc.density$x[which.max(ssc.density$y)]
    d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
    o<-order(d)
    sel.idx<-o[1:ncells]
    sel.xx<-xx[sel.idx,]
    cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), i, "\n");flush.console()
    plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", xlim=xlimi, bty="n", main="")
    text(0, mean(density(xx[,3])$y), paste(data.info[i, "effector"], data.info[i, "effector.level"], sep=" "), pos=4)
  }
}


wf<-flowSet(x)


sel.idx<-which(data.info[,"reg.prom"]=="HCE" & is.na(data.info[,"cofactor"]) & is.na(data.info[,"Enzyme"]) & is.na(data.info[,"substrate"]) )
sel.data.info<-data.info[sel.idx,]
myfplot(wf[sel.idx], sel.data.info, 3)


##X11()







##========= various nitrotoluene derivatives
##x<-list()
##cnt<-1
data.dir<-paste(root, "111107", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info2.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
##write.table(file.names, file="tmp.txt", quote=FALSE, sep="\n\r", row.names=F, col.names=F)
## -------------------------------------------
##x11()
##myfplot(wf, data.info, 3)


## Controls
sel.idx<-which(data.info[,"effector"]=="Control")
##sel.idx<-c(1:20)
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
myfplot(sel.data, sel.data.info, 3)



##========= Sensitivity comparison of Toluene, Benzene, Isoprene 
##x<-list()
##cnt<-1
data.dir<-paste(root, "111123", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
##write.table(file.names, file="tmp.txt", quote=FALSE, sep="\n\r", row.names=F, col.names=F)
## -------------------------------------------

## Controls
sel.idx<-which(data.info[,"effector"]=="Control" | data.info[,"effector"]=="Benzene")
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
elevel<-unlist(lapply(strsplit(sel.data.info[,"effector.level"], " "), function(x){x[1]}))
elevel[1]<-0
o<-order(elevel)
##X11()
##xx<-exprs(sel.data[[7]])
##plot(xx[,"FL1-H"])
myfplot(sel.data[o], sel.data.info[o,], 1)






##========= Sensitivity comparison of Isoprene 
##x<-list()
##cnt<-1
data.dir<-paste(root, "111123", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
##write.table(file.names, file="tmp.txt", quote=FALSE, sep="\n\r", row.names=F, col.names=F)
## -------------------------------------------

## Controls
sel.idx<-which(data.info[,"effector"]=="Control" | data.info[,"effector"]=="Isoprene")
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
elevel<-unlist(lapply(strsplit(sel.data.info[,"effector.level"], " "), function(x){x[1]}))
elevel[1]<-0
o<-order(elevel)
X11()
##plot(xx[,"FL1-H"])
myfplot(sel.data[o], sel.data.info[o,], 1)




## ============== 120228 various cell strain
##x<-list()
##cnt<-1
data.dir<-paste(root, "120228 various cell strain", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
## 
sel.idx<-c(9, 10, 6,7)
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
myfplot(sel.data, sel.data.info, 1, xlimi=c(0,800))

xx<-exprs(sel.data[[1]])



## ============== 121204 IspGESS promoter1
##x<-list()
##cnt<-1
data.dir<-paste(root, "121203 IspGESS promoter1", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
## 
sel.idx<-which(data.info[,"reg.prom"]=="TRC")
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
myfplot(sel.data, sel.data.info, 1, xlimi=c(0,800))


## ============== 121204 IspGESS promoter2
##x<-list()
##cnt<-1
data.dir<-paste(root, "121204 IspGESS promoter2", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
## 
sel.idx<-which(data.info[,"reg.prom"]=="ACE")
sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
myfplot(sel.data, sel.data.info, 1, xlimi=c(0,800))


## ============== 130122 B.subtilis candidates
##x<-list()
##cnt<-1
data.dir<-paste(root, "130122 B.subtilis candidates", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-rbind(data.info, tmp.info)
##data.info<-tmp.info
##wf<-flowSet(x)
## 
##sel.idx<-which(data.info[,"reg.prom"]=="ACE")
sel.data<-wf##[sel.idx]
sel.data.info<-data.info##[sel.idx,]
myfplot(sel.data, sel.data.info, 1, xlimi=c(0,800))






## ============== 20130307 Bindu sensitivity
x<-list()
cnt<-1
data.dir<-paste(root, "20130307", sep="")
file.names<-dir(data.dir, pattern="*.[0-9]$")
tmp.info<-as.matrix(read.csv(paste(data.dir, "/info.csv", sep="")))
for(i in 1:nrow(tmp.info)){
  x[[cnt]] <- read.FCS(paste(data.dir, "/", tmp.info[i,23], sep=""), transformation=FALSE)  
  cnt<-cnt+1 
}
data.info<-tmp.info
wf<-flowSet(x)
## 
##sel.idx<-which(data.info[,"reg.prom"]=="ACE")
sel.data<-wf##[sel.idx]
sel.data.info<-data.info##[sel.idx,]
x11()
##sel.idx<-2
##sel.data<-wf[sel.idx]
sel.data.info<-data.info[sel.idx,]
##png("foo.png")
myfplot(sel.data, sel.data.info, 1, xlimi=c(0,800))
savePlot(filename="plot.jpg", type="jpeg");
dev.off()

xx<-exprs(wf[[2]])
plot(density(xx[,"FSC-H"]), log="x")
plot(density(xx[,"SSC-H"])$y, log="x", type="l")
plot(density(log(xx[,"SSC-H"]))$y, type="l")












































 -------------------------------------------
## total
## FSC boxplot
## -------------------------------------------

par(ps=22, mfrow=c(1,nrow(data.info)+1), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0), bty="n")
xx<-exprs(wf[[1]])
plot(xx[,"FSC-H"], ylim=c(0, 550), type='n')
for(i in 1:length(wf)){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=0)
  xx<-xx[-del.idx,]
  if(i==1){
    boxplot(xx[,"FSC-H"], ylim=c(0, 550))
  }else{
    boxplot(xx[,"FSC-H"], ylim=c(0, 550), yaxt='n')
  }
}

## SSC boxplot
par(ps=22, mfrow=c(1,nrow(data.info)+1), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0), bty="n")
xx<-exprs(wf[[1]])
plot(xx[,"SSC-H"], ylim=c(0, 1000), type='n')
for(i in 1:length(wf)){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=0)
  xx<-xx[-del.idx,]
  if(i==1){
    boxplot(xx[,"SSC-H"], ylim=c(0, 1000))
  }else{
    boxplot(xx[,"SSC-H"], ylim=c(0, 1000), yaxt='n')
  }
}

## -------------------------------------------
##     control  
## -------------------------------------------
ctl.idx<-which(data.info[,16]=="Control")[-2]
ctl.info<-data.info[ctl.idx,]
write.table(ctl.info, file="tmp.csv", quote=FALSE, sep=",", row.names=F)


## FSC
par(ps=32, mfrow=c(1,length(ctl.idx)+1), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0), bty="n")
xx<-exprs(wf[[6]])
plot(xx[,"FSC-H"], ylim=c(0, 550), type='n')
for(i in ctl.idx){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=0)
  xx<-xx[-del.idx,]
  if(i==6){
    boxplot(xx[,"FSC-H"], ylim=c(0, 550))
  }else{
    boxplot(xx[,"FSC-H"], ylim=c(0, 550), yaxt='n')
  }
}


## SSC
par(ps=32, mfrow=c(1,length(ctl.idx)+1), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0), bty="n")
xx<-exprs(wf[[6]])
plot(xx[,"SSC-H"], ylim=c(0, 1000), type='n')
for(i in ctl.idx){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=0)
  xx<-xx[-del.idx,]
  if(i==6){
    boxplot(xx[,"SSC-H"], ylim=c(0, 1000))
  }else{
    boxplot(xx[,"SSC-H"], ylim=c(0, 1000), yaxt='n')
  }
}


## ----- 
par(ps=22, mfrow=c(3,2), mai=c(1, 1, 0, 0))
for(i in ctl.idx[c(1:3)]){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=50)
  xx<-xx[-del.idx,]
  ##select 10000 cells from
  
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  sel.idx<-o[1:10000]
  sel.xx<-xx[sel.idx,]
  cat(fsc.max, ssc.max, "dim:", dim(sel.xx), i, "\n");flush.console()
  plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  plot(density(sel.xx[,"FL1-H"]), col="#00000022", ylab="FL1", main="Scatter plot (FL1)")
}

## -------------------------------------------




















par(ps=28, mfrow=c(length(ctl.idx),2))
for(i in ctl.idx){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=0)
  xx<-xx[-del.idx,]
  plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC")
  gate.fsc<-c(10,400)
  gate.ssc<-c(550,800)
  idx<-which(xx[,"FSC-H"]>=gate.fsc[1] & xx[,"FSC-H"]<=gate.fsc[2] & xx[,"SSC-H"]>=gate.ssc[1] & xx[,"SSC-H"]<=gate.ssc[2])
  points(xx[idx,"FSC-H"], xx[idx,"SSC-H"], pch=19, col="#0000ff22")
  abline(v=gate.fsc[1], col="blue", lwd=2, lty=2)
  abline(v=gate.fsc[2], col="blue", lwd=2, lty=2)
  abline(h=gate.ssc[1], col="blue", lwd=2, lty=2)
  abline(h=gate.ssc[2], col="blue", lwd=2, lty=2)
  ##text(200, 1000, paste("FSC:", gate.fsc[1], "/", gate.fsc[2], sep=""))
  ##text(200, 900, paste("FSC:", gate.ssc[1], "/", gate.ssc[2], sep=""))
  length(idx)
  
  
  
  xxx<-density(xx[idx,"FL1-H"])
  plot(xxx, main="", xlim=c(-150, 700))
  polygon(xxx, col="#0000ff22", border="blue", main="")
  ##text(10, 0.004, c(paste("Mean:", round(mean(xx[idx,"FL1-H"]),1), sep="")))
  ##text(0, 0.0045, c(paste("Median:", round(median(xx[idx,"FL1-H"]),1), sep="")))
  ##text(30, 0.0040, paste("SD:", round(sd(xx[idx,"FL1-H"]),1), sep=""))
  abline(v=mean(xx[idx,"FL1-H"]), col="red", lwd=4, lty=2)
  abline(v=median(xx[idx,"FL1-H"]), col="blue", lwd=4, lty=2)
  
}






xx<-exprs(wf[[6]])

plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC")
gate.fsc<-c(500,700)
gate.ssc<-c(420,600)
idx<-which(xx[,"FSC-H"]>=gate.fsc[1] & xx[,"FSC-H"]<=gate.fsc[2] & xx[,"SSC-H"]>=gate.ssc[1] & xx[,"SSC-H"]<=gate.ssc[2])
points(xx[idx,"FSC-H"], xx[idx,"SSC-H"], pch=19, col="#0000ff22")
abline(v=gate.fsc[1], col="blue", lwd=2, lty=2)
abline(v=gate.fsc[2], col="blue", lwd=2, lty=2)
abline(h=gate.ssc[1], col="blue", lwd=2, lty=2)
abline(h=gate.ssc[2], col="blue", lwd=2, lty=2)
##text(200, 1000, paste("FSC:", gate.fsc[1], "/", gate.fsc[2], sep=""))
##text(200, 900, paste("FSC:", gate.ssc[1], "/", gate.ssc[2], sep=""))
length(idx)

xxx<-density(xx[idx,"FL1-H"])
plot(xxx, main="", xlim=c(-150, 700))
polygon(xxx, col="#0000ff22", border="blue", main="")
text(10, 0.004, c(paste("Mean:", round(mean(xx[idx,"FL1-H"]),1), sep="")))
text(0, 0.0045, c(paste("Median:", round(median(xx[idx,"FL1-H"]),1), sep="")))
text(30, 0.0040, paste("SD:", round(sd(xx[idx,"FL1-H"]),1), sep=""))
abline(v=mean(xx[idx,"FL1-H"]), col="red", lwd=4, lty=2)
abline(v=median(xx[idx,"FL1-H"]), col="blue", lwd=4, lty=2)
##}




###===========================================================



x<-list()
for(i in 1:length(file.names)){
  x[[i]] <- read.FCS(paste(data.dir, "/", file.names[i], sep=""), transformation=FALSE)  
}
wf<-flowSet(x)

##boxplot
par(ps=28, mfrow=c(1,nrow(data.info)), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0))
for(i in 1:length(file.names)){
  xx<-exprs(wf[[i]])
  if(i==2){
    boxplot(xx[,"FSC-H"], ylim=c(0, 850))
  }else{
    boxplot(xx[,"FSC-H"], ylim=c(0, 850), yaxt='n')
  }
}



###===========================================================



data.dir<-paste(root, "110830", sep="")
##file.names<-dir(data.dir, pattern="[0-9]$")
file.names<-dir(data.dir, pattern="GESS")
x<-list()
for(i in 1:length(file.names)){
  x[[i]] <- read.FCS(paste(data.dir, "/", file.names[i], sep=""), transformation=FALSE)  
}
wf<-flowSet(x)

##boxplot
par(ps=28, mfrow=c(1,nrow(data.info)), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0))
for(i in 1:length(file.names)){
  xx<-exprs(wf[[i]])
  if(i==2){
    boxplot(xx[,"FSC-H"], ylim=c(0, 850))
  }else{
    boxplot(xx[,"FSC-H"], ylim=c(0, 850), yaxt='n')
  }
}



par(ps=28, mfrow=c(2,nrow(data.info)), mai=c(0, 0, 0, 0), omi=c(0, 0, 0, 0))
for(i in 1:length(file.names)){

xx<-exprs(wf[[i]])
  
plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC")
gate.fsc<-c(10,400)
gate.ssc<-c(550,800)
idx<-which(xx[,"FSC-H"]>=gate.fsc[1] & xx[,"FSC-H"]<=gate.fsc[2] & xx[,"SSC-H"]>=gate.ssc[1] & xx[,"SSC-H"]<=gate.ssc[2])
points(xx[idx,"FSC-H"], xx[idx,"SSC-H"], pch=19, col="#0000ff22")
abline(v=gate.fsc[1], col="blue", lwd=2, lty=2)
abline(v=gate.fsc[2], col="blue", lwd=2, lty=2)
abline(h=gate.ssc[1], col="blue", lwd=2, lty=2)
abline(h=gate.ssc[2], col="blue", lwd=2, lty=2)
##text(200, 1000, paste("FSC:", gate.fsc[1], "/", gate.fsc[2], sep=""))
##text(200, 900, paste("FSC:", gate.ssc[1], "/", gate.ssc[2], sep=""))
length(idx)



xxx<-density(xx[idx,"FL1-H"])
plot(xxx, main="", xlim=c(-150, 700))
polygon(xxx, col="#0000ff22", border="blue", main="")
##text(10, 0.004, c(paste("Mean:", round(mean(xx[idx,"FL1-H"]),1), sep="")))
##text(0, 0.0045, c(paste("Median:", round(median(xx[idx,"FL1-H"]),1), sep="")))
##text(30, 0.0040, paste("SD:", round(sd(xx[idx,"FL1-H"]),1), sep=""))
abline(v=mean(xx[idx,"FL1-H"]), col="red", lwd=4, lty=2)
abline(v=median(xx[idx,"FL1-H"]), col="blue", lwd=4, lty=2)

}


































adata<-data.frame(media=c(rep("2step", 8), rep("casa", 8)), 
                  strain=c(rep("Celedx16", 4), rep("pPRO", 4), rep("Celedx15", 4), rep("pPRO", 4)), 
                  substrate=rep(c("none", "pNP", "pNPG2", "pNPG3"),4))


mData <- data.frame(labelDescription=c("media", "strain", "substrate"))

annod<-AnnotatedDataFrame()
##pData(wf)$name<-adata$media
##pData(wf)$strain<-adata$strain
##pData(wf)$substrate<-adata$substrate
##varLabels(wf)[1]<-"media"


par(ps=28, mfrow=c(1,2))
##for(i in 1:4){
  xx<-exprs(wf[[4]])
  plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC")
  gate.fsc<-c(500,700)
  gate.ssc<-c(420,600)
  idx<-which(xx[,"FSC-H"]>=gate.fsc[1] & xx[,"FSC-H"]<=gate.fsc[2] & xx[,"SSC-H"]>=gate.ssc[1] & xx[,"SSC-H"]<=gate.ssc[2])
  points(xx[idx,"FSC-H"], xx[idx,"SSC-H"], pch=19, col="#0000ff22")
  abline(v=gate.fsc[1], col="blue", lwd=2, lty=2)
  abline(v=gate.fsc[2], col="blue", lwd=2, lty=2)
  abline(h=gate.ssc[1], col="blue", lwd=2, lty=2)
  abline(h=gate.ssc[2], col="blue", lwd=2, lty=2)
  ##text(200, 1000, paste("FSC:", gate.fsc[1], "/", gate.fsc[2], sep=""))
  ##text(200, 900, paste("FSC:", gate.ssc[1], "/", gate.ssc[2], sep=""))
  length(idx)
  
  xxx<-density(xx[idx,"FL1-H"])
  plot(xxx, main="", xlim=c(-150, 700))
  polygon(xxx, col="#0000ff22", border="blue", main="")
  text(10, 0.004, c(paste("Mean:", round(mean(xx[idx,"FL1-H"]),1), sep="")))
  text(0, 0.0045, c(paste("Median:", round(median(xx[idx,"FL1-H"]),1), sep="")))
  text(30, 0.0040, paste("SD:", round(sd(xx[idx,"FL1-H"]),1), sep=""))
  abline(v=mean(xx[idx,"FL1-H"]), col="red", lwd=4, lty=2)
  abline(v=median(xx[idx,"FL1-H"]), col="blue", lwd=4, lty=2)
##}












hist(xx, freq = T, col = "grey")



## population
n=1000000
yy<-rnorm(n, mean=400, sd=1)
##yyy<-density(yy)
hist(yy, breaks=500, xlab="FL1-H", ylab="Count", main="")
##plot(yyy, xlab="FL1-H", ylab="Count", main="")
##hist

## sample
xx<-sample(yy, 100)
hist(xx, breaks=20, xlab="FL1-H", ylab="Count", main="", ps=20)














wf <- workFlow(ITN)
asinh <- arcsinhTransform()
tl <- transformList(colnames(ITN)[3:7], asinh, transformationId = "asinh")
add(wf, tl)


lg <- lymphGate(Data(wf[["asinh"]]), channels=c("SSC", "CD3"), preselection="CD4", filterId="TCells", eval=FALSE, scale=2.5)
add(wf, lg$n2gate, parent="asinh")
print(xyplot(SSC ~ CD3| PatientID, wf[["TCells+"]], par.settings=list(gate=list(col="red", fill="red", alpha=0.3))))





